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Overview: 


A  novel  approach  to  solving  large  multi-scale  image  based  problems  was  explored  opening  up  the 
possibility  of  solving  for  very  large  finite  element  and  finite  volume  problems  on  modest 
computational  platforms.  This  is  being  implemented  as  a  software  tool  as  part  of  'Material  Toolbox' 
for  use  by  the  research  and  development  community. 

Background: 

Synthetic  and  natural  micro-architectures  occur  frequently  (e.g.  MMC  foams,  bone,  etc.)  and 
multiphase  functionally  graded  composites  are  becoming  increasingly  popular  for  applications 
requiring  optimized/tailored  material  properties.  When  dealing  with  such  materials  computationally 
one  issue  which  immediately  arises  is  the  analysis  of  the  mechanical  properties  of  macroscopically 
inhomogeneous  multi-scale  structures.  The  bulk  response  of  these  structures  can  be  determined  by 
performing  'full'  finite  element  analysis,  that  is,  with  the  entire  geometry  discretized  at  a  resolution 
high  enough  to  accurately  model  the  smallest  length  scale  of  interest.  However,  these  full  models 
may  easily  exceed  hundreds  of  millions,  potentially  billions,  of  degrees  of  freedom  and  solving 
problems  of  this  magnitude  may  only  be  possible  with  the  use  of  supercomputing  facilities. 
Additionally,  in  an  iterative  optimization  process  where  the  'performance'  of  the  structure  may  be 
evaluated  thousands  of  times  the  use  of  full  FEA  simulations  becomes  highly  impractical.  When  the 
performance  of  a  structure  is  evaluated  in  an  optimization  process  typically  only  some  aspect  of  the 
bulk  response,  such  as  deflection,  is  considered.  For  such  properties  full  FEA  simulations  model  the 
problem  may  be  excessive.  In  the  present  project  a  novel  two  stage  approach  to  solving  such  large 
problem  by  performing  element  by  element  homogenization  of  the  micro-structure  followed  by 
solving  the  global  problem  with  a  coarser  mesh  as  will  be  outlined  below  was  explored. 

Approach  explored: 

The  approach  taken  is  effectively  based  on  creating  a  coarse  tetrahedral/hexahedral  discretization  of 
the  domain  using  traditional  volume  meshing  techniques  and  assigning  appropriate  material 
properties  based  on  a  finite  element  homogenization  based  on  high  resolution  mesh  at  the  micro- 
structural  level  of  the  macro  tetrahedra  or  hexahedra.  In  effect  two  length  scales  are  decoupled  by 
computing  effective  properties  using  the  finite  element  approach  for  each  macro-element.  The 
novelty  here  lies  in  effectively  discretizing  the  'full'  3D  mesh  into  larger  tetrahedra  and  hexahedra 
and  computing  homogenized  properties  for  each  macro-element  based  on  exact  meshed  domains 
representing  the  full  microstructural  complexity  within  the  macro-elements. 

Methodology/Proiect  Outline: 

The  project  was  implemented  in  several  steps: 

1.  Homogenization  for  linear  elastic  material  properties  on  a  hexahedral  domain 

In  the  limit  of  small  deformations,  elastic  materials  are  described  by  the  linear  theory  of  elasticity.  In 
an  elastic  body  at  equilibrium  and  in  the  absence  of  body  forces,  the  stress  and  strain  fields  a  and  e 
satisfy  the  equations 
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(1) 


1 

V  •  a  =  0,  e  =  -[Vu+(Vu)T], 

where  u(x)  :=  x  -  x0  is  the  displacement  of  the  point  x  of  the  body  relative  to  its  position  ;t0  in  a 
reference  stress-free  state  [Milton,  p.  22].  The  rank-two  symmetric  tensors  a  and  e  are  linked  by  the 
constitutive  relation 


a  =  C  e,  (2) 

where  C  is  the  (rank-four)  constitutive  tensor  (also  called  the  stiffness  tensor;  here  we  avoid  this 
term  due  to  possible  confusion  with  the  notion  of  a  stiffness  matrix  used  in  finite  element  analysis). 
In  a  heterogeneous  medium  C  is  position-dependent  [C  =  C(x)].  Thanks  to  this  relation,  eqs.  (1) 
can  be  reduced  to  a  single  (vectorial)  equation  for  the  displacement  u: 

V  •  {C  [Vu  +  (Vu)T]}  =  0.  (3) 

In  composite  media  one  can  usually  distinguish  three  distinct  length  scales  [Milton,  p.  7].  The 
microscale  is  characterized  by  the  typical  size  l±  of  heterogeneities:  it  is  the  scale  on  which  the 
intrinsic  material  properties,  such  as  C(x)  for  elasticity,  vary.  The  mesoscale  is  characterized  by  a 
length  l2  at  which  the  material  appears  "statistically  homogeneous".  At  this  scale  constitutive 
relations,  such  as  eq.  (2),  become  valid  in  the  sense  of  averages: 

a(x)  =  CeffO)  e(x),  (4a) 

where  Ceff(x)  is  the  effective  constitutive  tensor  at  point  x  and  the  overbar  denotes  the  moving 
average: 


a(x)  =  J  H(y  -  x)  e(y),  (4b) 

where  H(x)  is  some  convenient  window  function,  e.g.  H(x)  =  1  if  \x\  <  l2  and  0  otherwise. 

Finally,  the  macroscale  is  characterized  by  a  length  l3  at  which  the  effective  material  properties,  such 
as  Ceff(x),  undergo  substantial  variations.  If  the  mesoscale  is  much  smaller  than  the  macroscale 
(l2  «  l3),  it  becomes  possible  to  radically  simplify  the  conceptual  and  numerical  analysis  of  the 
composite  medium,  while  preserving  its  macroscopic  behaviour,  by  replacing  the  rapidly  varying 
constitutive  tensor  C(x)  by  the  much  more  slowly  varying  (or  constant)  tensor  Ceff(x). 

Equation  (4)  can  serve  as  the  operational  definition  of  the  effective  constitutive  tensor  Ceff  [Hashin 
1983].  To  calculate  its  value  for  a  given  composite  material,  one  needs  then  to  substitute  for  a(x) 
and  e(x)  valid  solutions  of  the  elasticity  equations  (1)— (2)  in  a  domain  H  being  the  support  of  the 
window  function  H(x).  The  tensor  Ceff  has  81  components,  only  21  of  which,  however,  are 
independent  [Saad  2009,  p.  80].  Owing  to  the  symmetry  properties  of  a,  e  and  Ceff  it  is  necessary  to 
use  six  linearly  independent  solutions  (ai(x),  (i  =  1;  2, ...,  6)  to  determine  all  components  of 

Ceff.  In  the  standard  shorthand  Voigt  notation,  in  which  the  stress  and  strain  tensors  are  written  as 
6-component  column  vectors  and  the  constitutive  tensor  as  a  6-by-6  symmetric  matrix,  the  final 
formula  for  Ceff  is 
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Cgff  =  [°i  ff2  ff6][ei  f2  e3  e4  c5  f6]_1,  (5) 

where  at  and  et  denote  column  vectors. 

Except  for  composites  with  very  particular  geometric  structure  (for  example  random  ensembles  of 
spherical  inclusions),  the  particular  solutions  (Oi(x),  e^x))  needed  to  calculate  the  effective 
constitutive  tensor  must  be  obtained  numerically,  and  this  is  the  approach  we  choose.  We  solve  eq. 
(4)  with  the  finite  element  method  in  a  hexahedral  domain  H.  We  use  ScanIP  to  mesh  the  domain 
with  linear  tetrahedral  elements  and  we  expand  each  component  of  the  unknown  field  u(x)  into 
globally  continuous,  piecewise  linear  scalar  basis  functions. 

The  choice  of  boundary  conditions  on  the  surface  dfl  of  H  is  an  important  aspect  of  numerical 
homogenization.  In  theory,  for  a  sufficiently  large  (i.e.  statistically  homogeneous)  domain  H  the 
same  value  of  the  effective  constitutive  tensor  Ceff  should  be  obtained  for  any  field  u(x)  satisfying 
eq.  (3)  in  fl,  and  hence  for  any  boundary  conditions  imposed  on  dfl.  In  practice,  however,  due  to 
limited  computational  resources,  simulations  are  often  done  on  domains  too  small  to  be  fully 
statistically  homogeneous.  In  this  case,  the  imposed  boundary  conditions  do  have  (some)  influence 
on  the  calculated  values.  Four  types  of  boundary  conditions  are  commonly  used  in  literature  [Ostoja- 
Starzewski  2006]: 

•  kinematic  uniform  boundary  conditions: 
u(x )  =  e0  •  x  for  x  E  dH, 

where  e0  is  a  constant  symmetric  tensor;  it  can  be  shown  that  e0  is  then  the  average  strain 
in  the  domain  H; 

•  static  uniform  boundary  conditions: 
a(X)  •  n(x)  =  ctq  •  n(x)  for  x  E  dH, 

where  ct0  is  a  constant  symmetric  tensor  and  n(x)  is  the  unit  vector  normal  to  the  surface 
dfi  at  x;  it  can  be  shown  that  a0  is  then  the  average  stress  in  the  domain  H; 

•  periodic  boundary  conditions: 

u(x  +  at )  =  u(x )  +  e0  •  x,  a(x  +  a*)  •  n(x  +  a*)  =  -d(x)  •  n(x )  for  x  E  dfL, 
for  at  being  one  of  three  linearly  independent  lattice  vectors  {a1;  a2,  a3}; 

•  mixed  orthogonal  boundary  conditions  (for  example,  uniform  strain  prescribed  on  two 
opposite  faces  of  a  cuboid  and  uniform  stress  prescribed  on  the  remaining  four  faces). 

In  each  case,  simulations  must  be  done  for  six  linearly  independent  values  of  the  constant  tensors  e0 
and  a0. 

On  the  boundary  dn,  the  uniform  boundary  conditions  introduce  artificial  features  that  would  be 
absent  in  a  domain  H  situated  within  a  larger  body  of  the  composite  material.  In  consequence, 
simulations  performed  using  kinematic  uniform  boundary  conditions  overestimate  the  effective 
constitutive  tensor  [Ostoja-Starzewski  2006],  while  the  application  of  static  uniform  boundary 
conditions  leads  to  an  underestimated  constitutive  tensor.  These  effects  are  particularly  strong 
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when  the  boundary  of  the  simulation  domain  dfi  intersects  material  interfaces.  Periodic  boundary 
conditions  do  not  have  this  deficiency,  but  only  as  long  as  they  are  applied  to  genuinely  periodic 
composites.  Mixed  orthogonal  boundary  conditions  are  only  applicable  to  materials  of  orthotropic  or 
higher  symmetry  [Hazanov  and  Amieur  1995],  and  the  edges  of  the  hexahedral  computational 
domain  fl  must  be  oriented  parallel  to  the  symmetry  axes  of  the  material.  If  these  requirements  are 
satisfied,  however,  mixed  orthogonal  boundary  conditions  provide  accurate  predictions  of  the 
effective  constitutive  tensor  already  for  small  computational  domains. 

Since  in  general  we  cannot  assume  the  presence  of  any  particular  symmetries  in  the  materials  to  be 
homogenized,  we  have  implemented  the  static  and  kinematic  uniform  boundary  conditions,  which 
are  applicable  to  all  composites  and  to  all  shapes  of  domains  H. 

2.  Homogenization  for  linear  elastic  material  properties  on  a  tetrahedral  domain 

The  homogenization  of  linear  elastic  material  properties  with  a  regular  hexahedral  domain,  as 
described  in  the  previous  section,  is  a  well-studied  topic.  However,  the  class  of  problems  of  interest 
to  this  project  (multi-scale  structures  with  irregular  domains)  required  that  a  different  approach  to 
the  calculation  of  the  effective  elastic  properties  be  taken.  The  more  straightforward  case  of  a 
regular  hexahedral  domain  may  be  addressed  by  dividing  the  domain  into  regular  hexahedral 
elements,  each  of  which  can  be  assigned  properties  using  the  classical  homogenization  techniques 
described  previously. 

For  the  more  general  case,  this  work  exploited  the  use  of  robust  all-tetrahedral  volume  meshing  as  a 
method  for  dividing  an  irregular  domain  into  smaller  sub-volumes  for  homogenization.  As  each  sub¬ 
volume  conformed  to  its  parent  macro-element  a  method  for  calculating  their  effective  properties 
was  developed. 

At  the  highest  level  the  developed  homogenization  process  involves  treating  the  sub-volume  as  if  it 
were  the  actual  macro-element.  Appropriate  boundary  conditions  are  applied  based  on  the  shape 
functions  of  the  macro-element  such  that  the  sub-volume  is  constrained  to  the  same  modes  of 
deformation.  A  series  of  finite  element  simulations  are  then  performed  in  order  to  determine  the 
sub-volume's  effective  properties.  A  more  detailed  description  of  this  process  in  2D  follows, 
however  the  approach  taken  extends  directly  to  3D  with  linear  tetrahedral  elements. 

3.  Constitutive  Matrix  Recovery  from  a  3  Noded  Triangle  Element 

To  highlight  the  basic  principle  of  the  developed  homogenization  method  this  section  will 
demonstrate  how  the  material  properties  of  a  simple  3-node  element  can  be  recovered  by  way  of 
virtual  testing.  This  in  itself  has  little  to  no  direct  practical  use,  but  the  principle  is  fundamental  to 
the  developed  homogenization  method. 

The  Constant  Strain  Triangle  (CST)  element  is  chosen  due  to  its  simplicity. 
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(*^m  ?  Urn  ) 


Figure  1:  A  3-node  constant  strain  triangle  element.  Nodes  are  numbered  anticlockwise. 

We  know  that  the  constitutive  matrix  and  element  geometry  determine  the  behavior  of  the 
element.  This  is  clear  from  how  the  element's  stiffness  matrix,  K,  is  calculated: 

K  =  tABTCB  (6) 


where  A  is  the  area  of  the  triangle,  C  the  constitutive  tensor  and  t  the  thickness  of  the  element.  The 
B  matrix  is  constructed  from  the  element's  shape  functions,  a  set  of  linear  displacement  functions. 
For  an  isotropic  material  in  2D  the  constitutive  tensor  is  defined  as: 


C  = 


1  —v2 


1  v 

V  1 
0  0 


0 

0 

1  -  v 
2  . 


(7) 


where  E  is  the  Young's  modulus  and  v  the  Poisson's  ratio.  Due  to  symmetry  in  the  matrix  there  are 
only  six  independent  components  and  hence  six  unknown  values  to  find  in  the  case  of  general 
anisotropy. 

We  know  that  applying  a  displacement  to  one  of  the  element's  nodes  will  result  in  a  force,  as 
described  by  Flooke's  Law: 


F  =  K  U  (8) 

where  F  is  the  force  vector  and  U  the  displacement  vector.  Substituting  equation  6  into  equation  8, 
it  is  possible  to  express  these  forces  in  terms  of  the  unknowns,  C\ 

F  =  tABTCBU  (9) 

The  expression  in  equation  9  will  provide  a  system  of  linear  equations  with  a  set  of  unknowns. 
Flowever,  for  a  single  test  the  system  is  underdetermined  and  it  can  be  shown  that  further  tests  are 
required  to  establish  an  over-determined  system.  In  2D  three  tests  are  required,  whereas  in  3D  six 
are  required.  A  solution  for  the  over-determined  system  is  then  calculated  using  the  least  squares 
method. 

Thus  far  it  has  been  shown  that  by  imposing  displacements  on  the  macro  element's  nodes  and 
measuring  the  resulting  forces  it  is  possible  to  recover  the  full  constitutive  matrix.  This  alone  is  of  no 
practical  use,  but  instead  provides  a  starting  point  for  bridging  the  gap  between  the  micro  and 
macro  scale. 
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In  the  previous  scheme  the  virtual  tests  could  be  described  using  a  simple  displacement  vector  as 
they  were  performed  on  a  single  element.  However,  in  this  case  the  virtual  tests  are  to  be  performed 
on  the  micro  mesh.  To  achieve  this  we  constrain  the  displacement  of  the  so-called  "external  micro 
nodes"  to  the  surface  of  the  macro  element.  This  constraint  is  used  so  that  the  area  discretized  by 
the  micro  elements  is  limited  to  behave  as  if  it  were  the  macro  element.  To  perform  the  virtual  tests 
we  first  compute  the  displacement  of  each  of  the  external  micro  nodes.  Given  a  displacement  vector 

U  =  [ut>  vt,  Uj,  Vj,  um,  vm } 

we  prescribe  the  external  micro  nodes  displacements  using  weightings  calculated  from  the  macro 
element's  shape  functions,  Nt,  Nj,  Nm: 


Ax  =  NiUi  +  NjUj  +  Nmum 
Ay  =  NiVt  +  NjVj  +  Nmvm 


Nodes  which  do  not  lie  on  the  macro  element's  boundary  are  left  unconstrained. 


Following  a  series  of  virtual  tests  a  number  of  forces  on  the  external  micro  nodes  will  have  been 
computed.  As  the  process  of  recovering  the  effective  constitutive  matrix  requires  the  forces  on  the 
macro  element  we  compute  the  effective  macro  forces  from  the  measured  micro  forces.  Similarly  to 
the  calculation  of  the  displacements,  the  effective  macro  forces  are  calculated  as  a  weighted  sum  of 
the  micro  forces  using  the  macro  element's  shape  function: 


Fi  = 


i  tk 


k= 1 


(11) 


where  Ft  is  the  macro  force  vector  for  macro  node  i,  n  the  number  of  external  micro  nodes  and  F£ 
and  F£  the  x  and  y  components  of  the  micro  force  vector  for  node  k.  Substituting  these  effective 
macro  forces  into  equation  9  we  are  able  to  compute  the  effective  constitutive  matrix  for  the  sub¬ 
volume. 


4.  Validation: 


As  part  of  the  validation  of  the  developed  homogenization  technique  for  tetrahedral  sub-volumes  a 
homogeneous  sub-volume  with  fully  anisotropic  material  properties  was  used  as  a  "sanity  test".  The 
input  material  properties  were  accurately  recovered. 

More  interestingly,  the  technique  was  also  used  to  recover  the  effective  properties  of  a  real  world 
structure  and  compared  to  the  results  obtained  using  classical  methods.  The  Schoen  Gyroid  was 
chosen  for  this  purpose  as  its  periodic  geometry  allows  periodic  boundary  conditions  to  be  used. 
These  boundary  conditions  are  often  considered  to  be  the  exact  solution. 

Figures  2  and  3  show  the  effective  Young's  modulus  and  shear  modulus  of  a  sample  of  the  Schoen 
Gyroid  structure  at  various  volume  fractions.  The  effective  properties  were  calculated  using  the 
classical  homogenization  technique  with  kinematic  uniform  boundary  conditions  ("KUBC"),  periodic 
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boundary  conditions  and  with  the  developed  method.  In  each  case  a  different  domain  was  used  to 
ensure  a  fair  comparison: 


•  KUBC- a  cube  of  8  x  8  x  8  unit  cells 

•  Periodic  -  a  single  unit  cell 

•  Tetrahedral  method  -  a  single  regular  tetrahedron  containing  83  unit  cells 


The  results  show  an  excellent  agreement  between  the  different  methods. 


Apparent  properties  -  8  x  8  x  8  unit  cells  (Schoen  Gyro  id) 


0.2  0  4  os  o.a 

Volume  fraction 


KUBC  — 
PERIODIC 

Tetrahedral  Method  + 
Voiqt  bound 
HS+ 


Figure  2:  Young's  modulus  of  the  Schoen  Gyroid  at  various  volume  fractions 
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Apparent  properties  -8x8x8  unit  cells  (Schoen  Gyroid) 
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Figure  3:  Shear  modulus  of  the  Schoen  Gyroid  at  various  volume  fractions 

5.  Sequential  Solution  of  Two  Stage  Problem 

Following  the  development  of  a  technique  to  determine  an  effective  constitutive  matrix  for  an 
arbitrary  tetrahedral  sub-volume  we  addressed  the  issue  of  multi-scale  problems.  Of  particular 


Page  |  7 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited 


interest  are  the  set  of  problems  having  an  irregular  (i.e.  non-cuboidal)  domain.  While  problems  of  a 
more  regular  nature  may  be  addressed  with  more  conventional  methods  of  determining  effective 
constitutive  matrices,  they  are  never  the  less  addressable  using  the  methods  developed  in  this  work. 

For  many  problems  it  is  highly  impractical  to  attempt  to  include  all  length  scales  in  a  finite  element 
model,  consequently  it  is  often  desirable  to  only  capture  the  coarser  details.  Rather  than  excluding 
the  smaller  length  scales  we  produce  a  coarse  mesh  with  appropriate  homogeneous  material 
properties.  The  previous  section  describes  the  process  of  computing  these  properties. 

The  process  of  generating  the  macro  mesh  is  outlined  in  Figure  4.  The  homogeneous  domain  should 
conform  to  the  bounds  of  the  original  domain  as  closely  as  possible,  representing  the  result  of  a 
'shrink  wrap1  operation. 


(a) 


(b) 


(c) 


Figure  4:  Femoral  head  modeling  (a)  input  dataset,  (b)coarse  tetrahedralisation  (c)  micro-element 

Each  of  the  macro  elements  in  the  generated  mesh  is  subsequently  homogenized  using  the 
developed  technique.  As  each  macro  element  is  considered  as  an  independent  sub-volume  the 
processing  may  occur  in  either  series  or  in  parallel,  depending  on  the  available  computational 
resources.  The  final  result  is  a  macroscopic  homogenous  model  with  varying  material  properties 
which  can  be  exported  to  a  traditional  finite  element  package. 

Conclusions: 

Project  showed  that  the  approach  could  be  successfully  used  on  problems  of  interest  for  the  elasto- 
statics  case  and  a  complete  pipeline  from  high  resolution  3D  image  data  through  to  multiscale  (two 
scale)  analysis.  The  approach  is  currently  being  implemented  in  code  for  release  in  a  version  of 
ScanIP.  This  is  being  carried  out  by  Dr  Wojciech  Smigaj  for  a  broad  range  of  physics  (rather  than 
solely  for  elasto-statics  problems  considered  under  the  auspices  of  EOARDS  grant). 

Dissemination: 

-Invited  abstract  accepted  for  oral  presentation  at  Telford-UKIERI  workshop  on 

"Anisotropic,  heterogeneous  and  cellular  materials:  From  microarchitecture  to  macro-level 
response"  in  Edinburgh,  December  12-14th,  2013.  See  attached  abstract. 

-Abstract  submitted  for  oral  presentation  at  2nd  International  Congress  on  3D  Materials  Science  2014 
on  "A  novel  approach  to  multiscale  homogenisation"  in  Annecy,  June  29th  -  July  2nd,  2014. 
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-Abstract  submitted  for  oral  presentation  at  11th  World  Congress  on  Computational  Mechanics 
(WCCM  XI)  on  "A  novel  approach  to  multiscale  homogenisation  for  3D  micro-structures"  in 
Barcelona,  July  20-25th,  2014. 

A-Paper  in  preparation  on  "Multi-scale  homogeneisation"  (  Formatted:  Not  Highlight 

-Visit  to  be  arranged  to  AFRL  Dayton,  Ohio  to  present  work  carried  out  and  explore  collaboration 
with  AFRL  research  staff 
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